/*************************************************************
** 	Incentive Effects of Recall Elections:  		        **
**	Evidence from Criminal Sentencing in California Courts 	**
** 	Analysis for submission.do								**
** 	Conduct all analysis for paper							**
** 	Latest revision: August 17, 2021				        **
************************************************************/

clear
set more off
set seed 10242004
set matsize 1000

use "recall_data_main.dta"

/********************************************************************************
*** 1. Main Effect 		   													  ***
********************************************************************************/
/********************************************************************************
*** 1a. Graphical Evidence (Figures 1 in paper and C.1 in SI)                 ***
********************************************************************************/
/********************************************************************************
*** 1a(1). Unadjusted             										      ***
********************************************************************************/
local vrange = "0.2(0.1)0.5"
local fitmodel = "lfit"
local fitmodel2 = "lpoly"
local dv = "sentence_normed_truncated"
local ar = 0.67
/* Petition Announced */
preserve
keep if abs(time_a)<=45
fastxtile n_all = time_a, nq(50)
gen n = .
replace n = n_all
collapse (mean) `dv' (median) time_a, by(n)
gen weight = max(0,(45-abs(time_a))/45)
twoway 	   (`fitmodel2' `dv' time_a if time_a <= 0 , lcolor(gs8) lpattern(solid)) ///
	       (`fitmodel2' `dv' time_a if time_a >=0 , lcolor(gs8) lpattern(solid)) ///
		   (`fitmodel' `dv' time_a if time_a <= 0 [aweight=weight], lcolor(maroon) lpattern(dash)) ///
	       (`fitmodel' `dv' time_a if time_a >=0 [aweight=weight], lcolor(maroon) lpattern(dash)) ///
		   (scatter `dv' time_a, mcolor(black) msymbol(circle_hollow)) ///
		    ,  name(graph_announce_unadjusted_paper, replace) ///
			title("Petition Announcement, Unadjusted")  xtitle("Date") ///
			xlabel( 0 "Jun 6" -45 "Apr 22" 45 "Jul 21",angle(45) nogrid) ///
			ytitle("Normalized Sentence") ylabel(`vrange', format(%3.1f)  nogrid) ///
			xline(0, lcolor(black) lpattern(dash)) legend(off) scheme(plotplain) aspectratio(`ar')
restore

/* Recall Election*/
preserve
keep if abs(time_r)<=45
fastxtile n_all = time_r, nq(50)
gen n = .
replace n = n_all
collapse (mean) `dv' (median) time_r, by(n)
gen weight = max(0,(45-abs(time_r))/45)
twoway 	   (`fitmodel2' `dv' time_r if time_r <= 0 , lcolor(gs8) lpattern(solid)) ///
	       (`fitmodel2' `dv' time_r if time_r >=0 , lcolor(gs8) lpattern(solid)) ///
		   (`fitmodel' `dv' time_r if time_r <= 0 [aweight=weight], lcolor(maroon) lpattern(dash)) ///
	       (`fitmodel' `dv' time_r if time_r >=0 [aweight=weight], lcolor(maroon) lpattern(dash)) ///
		   (scatter `dv' time_r, mcolor(black) msymbol(circle_hollow)) ///
		    ,  name(graph_recall_unadjusted_paper, replace) ///
			title("Recall Election, Unadjusted")  xtitle("Date") ///
			xlabel( 0 "Jun 5" -45 "Apr 21" 45 "Jul 20",angle(45)  nogrid) ///
			ytitle("Normalized Sentence") ylabel(`vrange', format(%3.1f) nogrid) ///
			xline(0, lcolor(black) lpattern(dash)) legend(off) scheme(plotplain) aspectratio(`ar')
restore 

/* For Appendix: Signature Certification*/
preserve
keep if abs(time_s)<=45
fastxtile n_all = time_s, nq(50)
gen n = .
replace n = n_all
collapse (mean) `dv' (median) time_s, by(n)
gen weight = max(0,(45-abs(time_s))/45)
twoway 	   (`fitmodel2' `dv' time_s if time_s <= 0 , lcolor(gs8) lpattern(solid)) ///
	       (`fitmodel2' `dv' time_s if time_s >=0 , lcolor(gs8) lpattern(solid)) ///
		   (`fitmodel' `dv' time_s if time_s <= 0 [aweight=weight], lcolor(maroon) lpattern(dash)) ///
	       (`fitmodel' `dv' time_s if time_s >=0 [aweight=weight], lcolor(maroon) lpattern(dash)) ///
		   (scatter `dv' time_s, mcolor(black) msymbol(circle_hollow)) ///
		    ,  name(graph_sigs_unadjusted_paper, replace) ///
			title("Signature Certification, Unadjusted")  xtitle("Date") ///
			xlabel(0 "Jan 24" -45 "Dec 10" 45 "March 10",angle(45)  nogrid) ///
			ytitle("Normalized Sentence") ylabel(`vrange', format(%3.1f) nogrid) ///
			xline(0, lcolor(black) lpattern(dash)) legend(off) scheme(plotplain) aspectratio(`ar')
restore 

/********************************************************************************
*** 1a(2). Residualized, using judge and statute fixed effects     			***
********************************************************************************/
/* Petition Announced */
preserve
drop if county == "Sacramento"
qui reg `dv' i.judge i.statute_code numberofcounts if abs(time_a)<=45,
qui predict ry if e(sample), resid
qui summ `dv' if e(sample), meanonly
qui replace ry = r(mean) + ry
qui reg time_a i.judge i.statute_code numberofcounts if abs(time_a)<=45
qui predict rx if e(sample), resid
qui summ time_a if e(sample), meanonly
qui replace rx = r(mean)+rx 

fastxtile n_all = rx, nq(50)
gen n = .
replace n = n_all

collapse (mean) ry (median) rx time_a, by(n)
gen weight = max(0,(45-abs(rx))/45)
twoway	   (`fitmodel2' ry rx if rx <= 0 , lcolor(gs8) lpattern(solid)) ///
	       (`fitmodel2' ry rx if rx  >=0 , lcolor(gs8) lpattern(solid)) ///
	       (`fitmodel' ry rx if rx <= 0 [aweight=weight], lcolor(maroon) lpattern(dash)) ///
	       (`fitmodel' ry rx if rx  >=0 [aweight=weight], lcolor(maroon) lpattern(dash)) ///
		   (scatter ry rx, mcolor(black) msymbol(circle_hollow)) ///
		    ,  name(graph_announce_adjusted_paper, replace) ///
			title("Petition Announcement, Adjusted")  xtitle("Date") ///
			xlabel( 0 "Jun 6" -45 "Apr 22" 45 "Jul 21",angle(45)  nogrid) ///
			ytitle("Normalized Sentence") ylabel(`vrange', format(%3.1f) nogrid) ///
			xline(0, lcolor(black) lpattern(dash)) legend(off) scheme(plotplain) aspectratio(`ar')			
restore

/* Recall Election*/
preserve
drop if county == "Sacramento"
qui reg `dv' i.judge i.statute_code numberofcounts if abs(time_r)<=45,
qui predict ry if e(sample), resid
qui summ `dv' if e(sample), meanonly
qui replace ry = r(mean) + ry
qui reg time_r i.judge i.statute_code numberofcounts if abs(time_r)<=45
qui predict rx if e(sample), resid
qui summ time_r if e(sample), meanonly
qui replace rx = r(mean)+rx 
summ rx
fastxtile n_all = rx, nq(50)
gen n = .
replace n = n_all

collapse (mean) ry (median) rx time_r, by(n)
gen weight = max(0,(45-abs(rx))/45)
twoway 	   (`fitmodel2' ry rx if rx <= 0 , lcolor(gs8) lpattern(solid)) ///
	       (`fitmodel2' ry rx if rx  >=0 , lcolor(gs8) lpattern(solid)) ///
	       (`fitmodel' ry rx if rx <= 0 [aweight=weight], lcolor(maroon) lpattern(dash)) ///
	       (`fitmodel' ry rx if rx  >=0 [aweight=weight], lcolor(maroon) lpattern(dash)) ///
			(scatter ry rx, mcolor(black) msymbol(circle_hollow)) ///
		    ,  name(graph_recall_adjusted_paper, replace) ///
			title("Recall Election, Adjusted")  xtitle("Date") ///
			xlabel( 0 "Jun 5" -45 "Apr 21" 45 "Jul 20",angle(45) nogrid) ///
			ytitle("Normalized Sentence") ylabel(`vrange', format(%3.1f) nogrid) ///
			xline(0, lcolor(black) lpattern(dash)) legend(off) scheme(plotplain) aspectratio(`ar')
restore


/*Generate Figure 1 in the paper, and export as pdf */
gr combine graph_announce_unadjusted_paper graph_recall_unadjusted_paper ///
		   graph_announce_adjusted_paper graph_recall_adjusted_paper, ///
		   imargin(0 0 0 0) graphregion(margin(l=22 r=22)) name(all_rd_plots, replace)
graph export "combined_rd_plots.pdf", replace

/* For Appendix: Signature Certification*/
preserve
drop if county == "Sacramento"
qui reg `dv' i.judge i.statute_code numberofcounts if abs(time_s)<=45,
qui predict ry if e(sample), resid
qui summ `dv' if e(sample), meanonly
qui replace ry = r(mean) + ry
qui reg time_s i.judge i.statute_code numberofcounts if abs(time_s)<=45
qui predict rx if e(sample), resid
qui summ time_s if e(sample), meanonly
qui replace rx = r(mean)+rx 
summ rx
fastxtile n_all = rx, nq(50)
gen n = .
replace n = n_all

collapse (mean) ry (median) rx time_s, by(n)
gen weight = max(0,(45-abs(rx))/45)
twoway 	   (`fitmodel2' ry rx if rx <= 0 , lcolor(gs8) lpattern(solid)) ///
	       (`fitmodel2' ry rx if rx  >=0 , lcolor(gs8) lpattern(solid)) ///
	       (`fitmodel' ry rx if rx <= 0 [aweight=weight], lcolor(maroon) lpattern(dash)) ///
	       (`fitmodel' ry rx if rx  >=0 [aweight=weight], lcolor(maroon) lpattern(dash)) ///
			(scatter ry rx, mcolor(black) msymbol(circle_hollow)) ///
		    ,  name(graph_sigs_adjusted_paper, replace) ///
			title("Signature Certification, Adjusted")  xtitle("Date") ///
			xlabel( 0 "Jan 24" -45 "Dec 10" 45 "March 10",angle(45) nogrid) ///
			ytitle("Normalized Sentence") ylabel(`vrange', format(%3.1f) nogrid) ///
			xline(0, lcolor(black) lpattern(dash)) legend(off) scheme(plotplain) aspectratio(`ar')
restore


/*Generate Figure C.1 in the SI, and export as pdf*/
gr combine graph_sigs_unadjusted_paper graph_sigs_adjusted_paper, ///
imargin(0 0 0 0) graphregion(margin(l=22 r=22)) name(sig_cert_appendix, replace)

graph export "signatures_rd_plots.pdf", replace

/********************************************************************************
*** 1b. RD with Optimal bandwidth      										  ***
*** This version of the do file loops over different                          ***
*** outcome variables; and restricting/not restricting cases                  ***
*** This generates estimates and exports them to the RD_template excel file   ***
*** Tables in paper generated: Table 1 (main results)                         ***
*** 						   Table 2(a) (normed_to_atc)                     ***
***                            Table 2(b) (charge bargain)                    ***
***                            Table B.1 (robust_untruncated_sentenc)         ***
***							   Table B.2 (robust_actual_sentence)			  ***
***							   Table B.3 (robust_onecount)                    ***
********************************************************************************/
set more off
local dvlist = "sentence_normed_truncated sentence_normed_truncated sentence_normed_truncated sentence_days sentence_normed_to_statute sentence_normed_to_atc_truncated charge_bargain"
local excellist = "main_results robust_filedates robust_onecount robust_actual_sentence robust_untruncated_sentenc normed_to_atc charge_bargain"
local treatlist = "announce recall"
local runninglist = "time_a time_r"
local n : word count `dvlist'
forvalues i = 1/`n' {
	local dv : word `i' of `dvlist'
    local sheet : word `i' of `excellist'
	matrix rd_results = J(6,4,.)
	if `i' == 2 {
		local cond_a = "& efiledate < mdy(6,6,2016)"
		local cond_r = "& efiledate < mdy(6,5,2018)"
		local extra_x = ""
		} 
	else if `i' == 3 {
		local cond_a = "& numberofcounts == 1"
		local cond_r = "& numberofcounts == 1"
		local extra_x = ""
	}
	else if `i' == 4 {
		local extra_x = "max_possible_sentence"
	}
	else if `i' == 5 {
		local cond_a = "& efiledate < mdy(6,6,2016)"
		local cond_r = "& efiledate < mdy(6,5,2018)"
		local extra_x = ""
	}
	else if `i' == `n' {
		local cond_a = "& plea_guilty == 1 & efiledate < mdy(6,6,2016)"
		local cond_r = "& plea_guilty == 1 & efiledate < mdy(6,5,2018)"
		local extra_x = ""
	}
		else {
		local cond_a = ""
		local cond_r = ""
		local extra_x = ""
		}
	forvalues j = 1/2 {
		local running : word `j' of `runninglist'
		local treat : word `j' of `treatlist'
		if `j' == 1 {
			local cond_x = "`cond_a'"
		}
		else if `j' == 2 {
			local cond_x = "`cond_r'"
		}
		
		/********************************************************************************
		*** 1b(1). Unadjusted        									  			  ***
		********************************************************************************/
		rdrobust `dv' `running' if abs(`running')<180 `cond_x', kernel(tri) vce(cluster cluster_county_statute) all
		matrix b_temp = e(beta_p_l)
		matrix v_temp = e(V_rb_l)
		matrix rd_results[1,2*`j'-1]=e(tau_bc)
		matrix rd_results[2,2*`j'-1]=e(se_tau_rb)
		matrix rd_results[3,2*`j'-1]= b_temp[1,1]
		matrix rd_results[4,2*`j'-1] = sqrt(v_temp[1,1])
		matrix rd_results[5,2*`j'-1]=e(h_l)
		matrix rd_results[6,2*`j'-1]=e(N_h_l)+e(N_h_r)
	
		/********************************************************************************
		*** 1b(2). Adjusted, Excludes Sacramento						  			  ***
		********************************************************************************/
		rdrobust `dv' `running' if abs(`running')<180 `cond_x', kernel(tri) vce(cluster cluster_judge_statute) ///
		covs(numberofcounts s_* j_*) all
		matrix b_temp = e(beta_p_l)
		matrix v_temp = e(V_rb_l)
		matrix rd_results[1,2*`j']=e(tau_bc)
		matrix rd_results[2,2*`j']=e(se_tau_rb)
		matrix rd_results[3,2*`j']= b_temp[1,1]
		matrix rd_results[4,2*`j'] = sqrt(v_temp[1,1])		
		matrix rd_results[5,2*`j']=e(h_l)
		matrix rd_results[6,2*`j']=e(N_h_l)+e(N_h_r)
		
		/********************************************************************************
		*** Including a large number of fixed effects via rdrobust can be computationally
		intensive. The original version of the paper reported the following alternative
		regression specification, which employs clustered standard errors, a triangular 
		kernel and MSE-optimal bandwidths drived from the unadjusted rdrobust model	  ***
		********************************************************************************/
		/*
		rdbwselect `dv' `running' if abs(`running')<180 & county !="Sacramento" `cond_x', kernel(tri) vce(cluster cluster_judge_statute)
		scalar bw = e(h_mserd)
		capture drop weight
		gen weight = max(0,(bw-abs(`running'))/bw)
		reg `dv' `treat'##c.`running' i.judge i.statute_code numberofcounts `extra_x' if abs(`running')<180 `cond_x' [pweight=weight], cluster(cluster_judge_statute)
		lincom 1.`treat'
		matrix rd_results[1,2*`j']=r(estimate)
		matrix rd_results[2,2*`j']=r(se)
		matrix rd_results[5,2*`j']=bw
		matrix rd_results[6,2*`j']=e(N)
		*/
		}
		capture drop rd_results*
		svmat rd_results
		/*Export to Excel table template*/
		export excel rd_results1-rd_results4 using "RD table template.xlsx", ///
		sheet("`sheet'") cell(C17) sheetmodify firstrow(variables)
}

capture drop rd_results*	  

/********************************************************************************
*** 1c. Balance on Charge Indicators (Petition Only) 					  	  ***
*** This robustness test checks for balance in the daily counts of charges	  ***
*** This generates Appendix Figure B.1										  ***
********************************************************************************/

set more off
matrix results = J(382,7,.)

forval t = 1(1)382{
	capture drop x dv
	disp in yellow "Charge: `t'"
	*count number of charges per day
	qui bys time_a: egen x=count(s_`t') if s_`t'==1
	qui replace x = 0 if mi(x)
	by time_a: egen dv = max(x)
	qui sum max_possible if s_`t'==1, d
	if missing("`r(max)'")==1{
		local max_poss = .
	}
	else if missing("`r(max)'")==0{
		local max_poss = `r(max)'
	}
	
	*estimate RD
	preserve
	keep time_a dv announce
	qui duplicates drop _all, force
	qui rdbwselect dv time_a if abs(time_a)<180, kernel(tri)
	scalar bw = e(h_mserd)
	capture drop weight
	qui gen weight = max(0,(bw-abs(time_a))/bw)
	capture qui reg dv announce##c.time_a [pweight=weight]
	restore

	*if not missing observations
	if missing("`e(N)'")==0{
		qui lincom 1.announce 
		matrix results[`t',1]=r(estimate)
		scalar p = 2*ttail(r(df),abs(r(estimate)/r(se)))
		matrix results[`t',2]=p
		matrix results[`t',3]= r(estimate) - invt(r(df),.975)*r(se)
		matrix results[`t',4]= r(estimate) + invt(r(df),.975)*r(se)
		matrix results[`t',5]=`t'
		matrix results[`t',6] = bw
		matrix results[`t',7] = `max_poss'
	} 
	else if missing("`e(N)'")==1 {
		disp in red "No Observations"
	}
	
}


capture drop results*
svmat results	  

drop x

*add noise to the max possible sentence so that all the results are not stacked against each other
gen x = log(results7+rnormal(0,50))
twoway (rspike  results3 results4 x,lcolor(maroon) horizontal) ///
		(scatter x results1, mcolor(gs6) msymbol(circle_hollow) msize(0.5)), ///
       yscale() name(gr_pvalues_unadj, replace) ytitle("Charge Severity (log days)", size(medium)) ///
	   xline(0, lcolor(black) lpattern(dash)) ///	   
	   xtitle("Estimate and 95% CI", size(medium)) legend(off) ///
	   ylabel(6.8(.5)9, angle(0) nogrid) xlab(,nogrid) ///
	   title("Charge-FE RD Estimates and 95% Confidence Intervals, Unadjusted",size(medium))  graphregion(color(white)) 
graph export "charge_balance_count.pdf", replace

/********************************************************************************
*** 2. Heterogeneous Effects (Petition Only)								  ***
*** 	Disaggregates by type of crime and race of defendant                  ***
*** This generates estimates and exports them to the RD_template excel file   ***
*** Tables in paper generated: Table 3 (by_crime_type)                        ***
***                            Table 4(A,B,C) (by_defendant_race)             ***
***                            Table 5 (by_defendant_race)                    ***
*** Figure in paper generated: Figure 3
********************************************************************************/
set more off
matrix rd_results = J(6,6,.)
local dv = "sentence_normed_truncated"
local excellist = "by_crime_type by_defendant_race"
local crimelist = "sex_crime violent_nonsex_667 nonviolent_crime_667"
local racelist  = "Black Hispanic White"
local running = "time_a"
local treat = "announce"
forvalues i = 1/2 {
	local sheet : word `i' of `excellist'
	if `i' == 1 {
		local conditionlist = "`crimelist'"
	}
	else {
		local conditionlist = "`racelist'"
	}
	local n : word count `conditionlist'
	forvalues j = 1/`n' {
		local cond : word `j' of `conditionlist'
		local cond = "`cond'" + " == 1"
		/********************************************************************************
		*** 2a. Unadjusted        									  			  ***
		********************************************************************************/
		rdrobust `dv' `running' if abs(`running')<180 & `cond', ///
		kernel(tri) vce(cluster cluster_county_statute) all
		matrix b_temp = e(beta_p_l)
		matrix v_temp = e(V_cl_l)
		matrix rd_results[1,2*`j'-1]=e(tau_bc)
		matrix rd_results[2,2*`j'-1]=e(se_tau_rb)
		matrix rd_results[3,2*`j'-1]= b_temp[1,1]
		matrix rd_results[4,2*`j'-1] = sqrt(v_temp[1,1])
		matrix rd_results[5,2*`j'-1]=e(h_l)
		matrix rd_results[6,2*`j'-1]=e(N_h_l)+e(N_h_r)
	
		/********************************************************************************
		*** 2b. Adjusted, Excludes Sacramento						  			  ***
		********************************************************************************/
		rdrobust `dv' `running' if abs(`running')<180 & `cond', ///
		kernel(tri) vce(cluster cluster_judge_statute) all covs(numberofcounts s_* j_*)
		matrix b_temp = e(beta_p_l)
		matrix v_temp = e(V_cl_l)
		matrix rd_results[1,2*`j']=e(tau_bc)
		matrix rd_results[2,2*`j']=e(se_tau_rb)
		matrix rd_results[3,2*`j']= b_temp[1,1]
		matrix rd_results[4,2*`j'] = sqrt(v_temp[1,1])
		matrix rd_results[5,2*`j']=e(h_l)
		matrix rd_results[6,2*`j']=e(N_h_l)+e(N_h_r)
		/********************************************************************************
		*** Including a large number of fixed effects via rdrobust can be computationally 
		intensive. The original version of the paper reported the following alternative
		regression specification, which employs clustered standard errors, a triangular 
		kernel and MSE-optimal bandwidths drived from the unadjusted rdrobust model	  ***
		********************************************************************************/
		/*
		rdbwselect `dv' `running' if abs(`running')<180 & county !="Sacramento" & `cond', kernel(tri) vce(cluster cluster_judge_statute)
		scalar bw = e(h_mserd)
		capture drop weight
		gen weight = max(0,(bw-abs(`running'))/bw)
		reg `dv' `treat'##c.`running' i.judge i.statute_code numberofcounts `extra_x' if abs(`running')<180 `cond_x' [pweight=weight], cluster(cluster_judge_statute)
		lincom 1.`treat'
		matrix rd_results[1,2*`j']=r(estimate)
		matrix rd_results[2,2*`j']=r(se)
		matrix rd_results[5,2*`j']=bw
		matrix rd_results[6,2*`j']=e(N)
		*/
		}
		capture drop rd_results*
		svmat rd_results
		/*Export to Excel table template*/
		export excel rd_results1-rd_results6 using "RD table template.xlsx", ///
		sheet("`sheet'") cell(C17) sheetmodify firstrow(variables)
	}
	
		/***********************************************************************************
		*** 2c. Hypothesis Tests of race Coefficient Differences						 ***
		***********************************************************************************/
		matrix H_results = J(12,3,.)
		/*Unadjusted Estimates*/
		qui rdrobust sentence_normed_truncated time_a if abs(time_a)<180 & White == 1, ///
		kernel(tri) vce(cluster cluster_county_statute) all
		local white_mean =e(tau_bc)
		local white_se = (se_tau_rb)
		matrix bpl = e(beta_p_l)
		local white_lsi = bpl[1,1]
		matrix vl = e(V_rb_l)
		local white_vpl = vl[1,1]	
		qui rdrobust sentence_normed_truncated time_a if abs(time_a)<180 & Black == 1, ///
		kernel(tri) vce(cluster cluster_county_statute) all
		local black_mean =e(tau_bc)
		local black_se = (se_tau_rb)
		matrix bpl = e(beta_p_l)
		local black_lsi = bpl[1,1]
		matrix vl = e(V_rb_l)
		local black_vpl = vl[1,1]
		qui rdrobust sentence_normed_truncated time_a if abs(time_a)<180 & Hispanic == 1, ///
		kernel(tri) vce(cluster cluster_county_statute) all
		local hispanic_mean =e(tau_bc)
		local hispanic_se = (se_tau_rb)
		matrix bpl = e(beta_p_l)
		local hispanic_lsi = bpl[1,1]		
		matrix vl = e(V_rb_l)
		local hispanic_vpl = vl[1,1]
		matrix H_results[1,1]=abs(`white_mean'-`black_mean')
		matrix H_results[2,1]=sqrt((`white_se')^2+(`black_se')^2)
		matrix H_results[3,1]=abs(`hispanic_mean'-`black_mean')
		matrix H_results[4,1]=sqrt((`hispanic_se')^2+(`black_se')^2)
		matrix H_results[5,1]=abs(`white_mean'-`hispanic_mean')
		matrix H_results[6,1]=sqrt((`hispanic_se')^2+(`white_se')^2)
		matrix H_results[7,1]=abs(`white_lsi'-`black_lsi')
		matrix H_results[8,1]=sqrt(`white_vpl'+`black_vpl')
		matrix H_results[9,1]=abs(`black_lsi'-`hispanic_lsi')
		matrix H_results[10,1]= sqrt(`black_vpl'+`hispanic_vpl')
		matrix H_results[11,1]=abs(`white_lsi'-`hispanic_lsi')
		matrix H_results[12,1]=sqrt(`white_vpl'+`hispanic_vpl')

		/*Adjusted Estimates*/		
		/*Unadjusted Estimates*/
		qui rdrobust sentence_normed_truncated time_a if abs(time_a)<180 & White == 1, ///
		kernel(tri) vce(cluster cluster_judge_statute) all covs(numberofcounts s_* j_*)
		local white_mean =e(tau_bc)
		local white_se = (se_tau_rb)
		matrix bpl = e(beta_p_l)
		local white_lsi = bpl[1,1]
		matrix vl = e(V_rb_l)
		local white_vpl = vl[1,1]	
		qui rdrobust sentence_normed_truncated time_a if abs(time_a)<180 & Black == 1, ///
		kernel(tri) vce(cluster cluster_judge_statute) all covs(numberofcounts s_* j_*)
		local black_mean =e(tau_bc)
		local black_se = (se_tau_rb)
		matrix bpl = e(beta_p_l)
		local black_lsi = bpl[1,1]
		matrix vl = e(V_rb_l)
		local black_vpl = vl[1,1]
		qui rdrobust sentence_normed_truncated time_a if abs(time_a)<180 & Hispanic == 1, ///
		kernel(tri) vce(cluster cluster_judge_statute) all covs(numberofcounts s_* j_*)
		local hispanic_mean =e(tau_bc)
		local hispanic_se = (se_tau_rb)
		matrix bpl = e(beta_p_l)
		local hispanic_lsi = bpl[1,1]		
		matrix vl = e(V_rb_l)
		local hispanic_vpl = vl[1,1]
		matrix H_results[1,3]=abs(`white_mean'-`black_mean')
		matrix H_results[2,3]=sqrt((`white_se')^2+(`black_se')^2)
		matrix H_results[3,3]=abs(`hispanic_mean'-`black_mean')
		matrix H_results[4,3]=sqrt((`hispanic_se')^2+(`black_se')^2)
		matrix H_results[5,3]=abs(`white_mean'-`hispanic_mean')
		matrix H_results[6,3]=sqrt((`hispanic_se')^2+(`white_se')^2)
		matrix H_results[7,3]=abs(`white_lsi'-`black_lsi')
		matrix H_results[8,3]=sqrt(`white_vpl'+`black_vpl')
		matrix H_results[9,3]=abs(`black_lsi'-`hispanic_lsi')
		matrix H_results[10,3]= sqrt(`black_vpl'+`hispanic_vpl')
		matrix H_results[11,3]=abs(`white_lsi'-`hispanic_lsi')
		matrix H_results[12,3]=sqrt(`white_vpl'+`hispanic_vpl')		
		svmat H_results
		export excel H_results1-H_results3 using "RD table template.xlsx", ///
		sheet("by_defendant_race") cell(T16) sheetmodify firstrow(variables)
		capture drop H_results*
		/************************************************************************************
		*** 2d. Racial disparities in charge severity, pre, post, total (Table 5 in paper)***
		************************************************************************************/
		eststo clear
		eststo m1: reg max_possible_sentence Black Hispanic  if time_a < 0  & (White==1|Black==1|Hispanic==1), cluster(county)
		eststo m2: xtreg max_possible_sentence Black Hispanic if time_a < 0  & (White==1|Black==1|Hispanic==1), fe i(judge) cluster(judge)
		eststo m3: reg max_possible_sentence Black Hispanic  if time_a >= 0  & (White==1|Black==1|Hispanic==1), cluster(county)
		eststo m4: xtreg max_possible_sentence Black Hispanic if time_a >= 0  & (White==1|Black==1|Hispanic==1), fe i(judge) cluster(judge)
		eststo m5: reg max_possible_sentence Black Hispanic if (White==1|Black==1|Hispanic==1), cluster(county)
		eststo m6: xtreg max_possible_sentence Black Hispanic if (White==1|Black==1|Hispanic==1) , fe i(judge) cluster(judge)
		estout * using "max_sentence_by_race.csv", cells(b(fmt(1)) se(par fmt(1))) stats(N) style(tex) replace
		/********************************************************************************
		*** 2e. Racial disparities in total sentencing	(Figure 3 in paper)           ***
		********************************************************************************/
		matrix results = J(6,4,.)
		local condlist = "White Black Hispanic"
		forval i = 1(1)3 {
			local cond: word `i' of `condlist'
			rdrobust log10_sentence_days time_a if abs(time_a)<180 & `cond' == 1, ///
				     kernel(tri) vce(cluster cluster_county_statute) all
			matrix vl = e(V_rb_l)
			matrix vr = e(V_rb_r)
			matrix results[2*`i'-1,1] =  e(tau_bc_l)
			matrix results[2*`i'-1,2] =  e(tau_bc_l) - 1.959964*sqrt(vl[1,1])
			matrix results[2*`i'-1,3] =  e(tau_bc_l) + 1.959964*sqrt(vl[1,1])
			matrix results[2*`i'-1,4]=0
			matrix results[2*`i',1] =  e(tau_bc_r)
			matrix results[2*`i',2] = e(tau_bc_r) - 1.959964*sqrt(vr[1,1])
			matrix results[2*`i',3] = e(tau_bc_r) + 1.959964*sqrt(vr[1,1])
			matrix results[2*`i',4]=1
			}
		capture drop results*	
		svmat results
		capture drop r
		gen r = "White" in 1/2
		replace r = "Black" in 3/4
		replace r = "Hispanic" in 5/6
		twoway (rspike results2 results3 results4) (scatter results1 results4) (line results1 results4), ///
		       by(r, col(3) legend(off) note("")) xtitle("") ytitle("Sentence Length (Days)") ///
		       xlabel(0 "Pre" 1 "Post") ylabel(2 "100" 2.39794 "250" 2.69897 "500" 2.8750613 "750")
		graph export "sentence_length_by_race.pdf", replace 	
		
/********************************************************************************
*** 3. Additional Robustness Checks					         				  ***
********************************************************************************/
/********************************************************************************
*** 3a. Placebo Tests -- Petition Only (Figure 2 in paper)                    ***
********************************************************************************/
/********************************************************************************
*** 3a(1). No covariates 			      								      ***
********************************************************************************/
preserve
set more off
set matsize 2000
local startdate = mdy(1,1,2016)
local enddate = mdy(12,31,2016)
scalar length = 1+`enddate'-`startdate'
matrix results = J(length,6,.)

forval t = `startdate'(1)`enddate'{
	disp in yellow "Day: `t'"
	capture drop time_x dummy 
	gen time_x = edate - `t'
	gen dummy = edate >= `t'
	qui rdrobust sentence_normed_truncated time_x if abs(time_x)<180, ///
			  kernel(tri) vce(cluster cluster_county_statute) all
	local t0 = `t'-`startdate'+1
	matrix results[`t0',1]= e(tau_bc)
	matrix results[`t0',2]=2*normal(-abs(e(tau_bc)/e(se_tau_rb)))
	matrix results[`t0',3]= e(tau_bc) - 1.959964*e(se_tau_rb)
	matrix results[`t0',4]= e(tau_bc) + 1.959964*e(se_tau_rb)
	matrix results[`t0',5]=`t'
	matrix results[`t0',6] = e(b_l)
}

capture drop results*
svmat results	  
format results5 %tdMon_dd 
cumul results2, gen(cumul_unadjusted)

qui summ results2 if results5 == mdy(6,6,2016)
local arrow_x = r(mean)
twoway (line cumul_unadj cumul_unadj) ///
	   (scatter cumul_unadj results2 if results5 !=mdy(6,6,2016), mcolor(gs6) msymbol(circle_hollow)) ///
       (scatter cumul_unadj results2  if results5 == mdy(6,6,2016), mcolor(maroon) msymbol(diamond) msize(medium) ) ///
	   (pcarrowi 0.3 `arrow_x' 0.05 `arrow_x' (1), mcolor(maroon) lcolor(maroon) mfcolor(maroon) barbsize(1) msize(2)), text(0.34 0.025 "Actual", color(maroon)) ///
       yscale() name(gr_pvalues_unadj, replace) xtitle("Placebo p-value", size(small)) ///
	   ytitle("Empirical Cumulative Distribution", size(small)) legend(off) ///
	   ylab(,nogrid format(%3.1f)) xlab(,nogrid) ///
	   title("Placebo RD P-Values, Unadjusted",size(small))
capture drop n
sort cumul_unadj
qui count if results1 !=.
local n_temp = r(N)
gen n = _n
qui summ n if results5==mdy(6,6,2016)
local p_rank = r(mean)
disp (`n_temp'-`p_rank')/(`n_temp'-1)

local arrow_x =mdy(6,6,2016)
local limit = 45
twoway (rspike results3 results4 results5 if results5!=20611 & abs(results5-mdy(6,6,2016))<`limit',lcolor(gs11)) ///
		(scatter results1 results5 if results5!=20611 & abs(results5-mdy(6,6,2016))<`limit', mcolor(gs6) msymbol(circle_hollow)) ///
		(scatter results1 results5 if results5==20611 & abs(results5-mdy(6,6,2016))<`limit', mcolor(maroon) msymbol(diamond) ) ///
		(rspike results3 results4 results5 if results5==20611 & abs(results5-mdy(6,6,2016))<`limit',lcolor(maroon) lwidth(medthick)) ///
		(pcarrowi -0.1 `arrow_x' -0.02 `arrow_x' (1), mcolor(maroon) lcolor(maroon) mfcolor(maroon) barbsize(1) msize(2)), text(-.11 `arrow_x' "Actual", color(maroon)) ///
		xlabel(,angle(45) nogrid) yline(0,lcolor(black)) legend(off) name(gr_ts_unadj, replace) xtitle("Date",size(small)) ///
		ylab(,format(%3.1f) nogrid) ytitle("Estimate and 95% CI", size(small)) ///
		title("Placebo RD Estimates and 95% Confidence Intervals, Unadjusted",size(small)) 
restore

/*
*NOTE: Section 3b takes approximately 52 hours to run (tested on Intel(R) Core(TM) i5-4590S CPU processor with 16.0 GB RAM); 	
/********************************************************************************
*** 3a(2). With covariates 			      								      ***
********************************************************************************/

preserve
set more off
set matsize 2000
local startdate = mdy(1,1,2016)
local enddate = mdy(12,31,2016)
scalar length = 1+`enddate'-`startdate'
matrix results = J(length,6,.)

forval t = `startdate'(1)`enddate'{
	disp in yellow "Day: `t'"
	capture drop time_x dummy 
	gen time_x = edate - `t'
	gen dummy = edate >= `t'
	qui rdrobust sentence_normed_truncated time_x if abs(time_x)<180, ///
			  kernel(tri) vce(cluster cluster_county_statute) covs(numberofcounts s_* j_*) all
	local t0 = `t'-`startdate'+1
	matrix results[`t0',1]= e(tau_bc)
	matrix results[`t0',2]=2*normal(-abs(e(tau_bc)/e(se_tau_rb)))
	matrix results[`t0',3]= e(tau_bc) - 1.959964*e(se_tau_rb)
	matrix results[`t0',4]= e(tau_bc) + 1.959964*e(se_tau_rb)
	matrix results[`t0',5]=`t'
	matrix results[`t0',6] = e(b_l)
	
}
local ar = 0.67
capture drop results*
svmat results	  
format results5 %tdMon_dd 
capture drop cumul*
cumul results2, gen(cumul)
local ar = 0.67
qui summ results2 if results5 == mdy(6,6,2016)
local arrow_x = r(mean)
twoway (line cumul cumul) ///
	   (scatter cumul results2 if results5 !=mdy(6,6,2016), mcolor(gs6) msymbol(circle_hollow)) ///
       (scatter cumul results2  if results5 == mdy(6,6,2016), mcolor(maroon) msymbol(diamond) msize(medium) ) ///
	   (pcarrowi 0.3 `arrow_x' 0.05 `arrow_x' (1), mcolor(maroon) lcolor(maroon) mfcolor(red) barbsize(1) msize(2)), text(0.33 0.025 "Actual", color(maroon)) ///
       yscale() name(gr_pvalues_adj, replace) xtitle("Placebo p-value",size(medium)) ///
	   ytitle("Empirical Cumulative Distribution",size(medium)) legend(off) aspectratio(`ar') ///
	   ylab(,nogrid format(%3.1f)) xlab(,nogrid) ///
	   title("Placebo RD P-Values, Adjusted",size(medium))

capture drop results*
svmat results	  
format results5 %tdMon_dd 
local arrow_x = mdy(6,6,2016)
twoway (rspike results3 results4 results5 if results5!=20611,lcolor(gs11)) ///
		(scatter results1 results5 if results5!=20611, mcolor(gs6) msymbol(circle_hollow)) ///
		(scatter results1 results5 if results5==20611, mcolor(maroon) msymbol(diamond) ) ///
		(rspike results3 results4 results5 if results5==20611,lcolor(maroon) lwidth(medthick)) ///
		(pcarrowi 0.35 `arrow_x' 0.22 `arrow_x' (1), mcolor(maroon) lcolor(maroon) mfcolor(maroon) barbsize(1) msize(2)), text(0.38 `arrow_x' "Actual", color(maroon)) ///
		xlabel(,angle(45) nogrid) yline(0,lcolor(black)) legend(off) name(gr_ts_adj, replace) xtitle("Date",size(medium)) /// 
		ylab(-0.4(0.1)0.2,format(%3.1f) nogrid) ytitle("Estimate and 95% CI",size(medium)) aspectratio(`ar') ///
		title("Placebo RD Estimates and 95% Confidence Intervals, Adjusted",size(medium)) 
capture drop n
sort cumul
qui count if results1 !=.
local n_temp = r(N)
gen n = _n
qui summ n if results5==mdy(6,6,2016)
local p_rank = r(mean)
disp (`n_temp'-`p_rank')/(`n_temp'-1)		

gr combine gr_pvalues_unadj gr_ts_unadj gr_pvalues_adj gr_ts_adj, xsize(7.5) ysize(5)
graph export "placebo_tests_both.pdf", replace		
restore
*/

/********************************************************************************
*** 3b. Tests for Bandwidth Artifacts - Petition Only (Figure B.2 in SI)      ***
********************************************************************************/
matrix bw_results_unadj = J(81,4,.)
matrix bw_results_adj = J(81,4,.)
forval bw = 10(1)90 {
	rdrobust sentence_normed_truncated time_a if abs(time_a)<180 , h(`bw' `bw')  b(`bw' `bw') ///
			  kernel(tri) vce(cluster cluster_county_statute) all
	matrix Vl = e(V_rb_l)
	matrix Vr = e(V_rb_r)
	local se = sqrt(Vl[1,1]+Vr[1,1])
	matrix bw_results_unadj[`bw'-9,1]=`bw'
	matrix bw_results_unadj[`bw'-9,2]=e(tau_bc) 
	matrix bw_results_unadj[`bw'-9,3]=e(tau_bc) - 1.959964*`se'
	matrix bw_results_unadj[`bw'-9,4]=e(tau_bc) + 1.959964*`se'

	rdrobust sentence_normed_truncated time_a if abs(time_a)<180 , h(`bw' `bw')  b(`bw' `bw') ///
			  covs(numberofcounts s_* j_*) kernel(tri) vce(cluster cluster_county_statute) all
	matrix Vl = e(V_rb_l)
	matrix Vr = e(V_rb_r)
	local se = sqrt(Vl[1,1]+Vr[1,1])
	matrix bw_results_adj[`bw'-9,1]=`bw'
	matrix bw_results_adj[`bw'-9,2]=e(tau_bc) 
	matrix bw_results_adj[`bw'-9,3]=e(tau_bc) - 1.959964*`se'
	matrix bw_results_adj[`bw'-9,4]=e(tau_bc) + 1.959964*`se'
}	
capture drop bw_results*
svmat bw_results_unadj
svmat bw_results_adj
local ar = 0.67
twoway (line bw_results_unadj3 bw_results_unadj1, lcolor(black) lpattern(dot)) ///
	   (line bw_results_unadj4 bw_results_unadj1, lcolor(black) lpattern(dot)) ///
	   (line bw_results_unadj2 bw_results_unadj1, lcolor(black) lpattern(solid)), ///
	   legend(off) name(gr_band_unadj, replace) xtitle("Bandwidth") ytitle("Estimate and 95% CI") /// 
	   xlab(10(10)90, nogrid) ylab(-0.3(0.1)0.5, nogrid format(%3.1f)) yline(0,lpattern(solid)) xline(45.1,lpattern(solid)) ///
	   aspectratio(`ar') title("Unadjusted Estimates, Varying Bandwidth") 

twoway (line bw_results_adj3 bw_results_adj1, lcolor(black) lpattern(dot)) ///
	   (line bw_results_adj4 bw_results_adj1, lcolor(black) lpattern(dot)) ///
	   (line bw_results_adj2 bw_results_adj1, lcolor(black) lpattern(solid)), ///
	   legend(off) name(gr_band_adj, replace) xtitle("Bandwidth") ytitle("Estimate and 95% CI") /// 
	   xlab(10(10)90, nogrid) ylab(-0.3(0.1)0.5, nogrid format(%3.1f)) yline(0,lpattern(solid)) xline(44.6,lpattern(solid)) ///
	   aspectratio(`ar') title("Adjusted Estimates, Varying Bandwidth") 

gr combine gr_band_unadj gr_band_adj, imargin(0 0 0 0) graphregion(margin(l=22 r=22)) xsize(7) ysize(3.5)
graph export "bandwidth_tests.pdf", replace		



/********************************************************************************
*** 4. Calculate Aggregate Effects                                           ***
********************************************************************************/
capture drop agg_results*
matrix agg_results = J(5,3,.)
/********************************************************************************
*** 4a. Assume LATE is ATE                                                    ***
********************************************************************************/
/********************************************************************************
*** 4a(1). Unadjusted                                                         ***
********************************************************************************/
summ max_possible if announce == 1 & abs(time_a)<45
scalar N_window = r(N)
scalar mean_window = r(mean)
rdrobust sentence_days time_a if abs(time_a)<180, kernel(tri) vce(cluster cluster_county_statute) all
/*Note -- these results are on the normalized scale -- adjusted later in program*/
matrix agg_results[1,1] = (e(tau_bc)-invnorm(.975)*e(se_tau_rb))*N_window/365
matrix agg_results[1,2] = e(tau_bc)*N_window/365
matrix agg_results[1,3] = (e(tau_bc)+invnorm(.975)*e(se_tau_rb))*N_window/365

/********************************************************************************
*** 4a(2). Adjusted                                                           ***
********************************************************************************/
rdrobust sentence_days time_a if abs(time_a)<180, kernel(tri) /// 
		 vce(cluster cluster_judge_statute) covs(numberofcounts s_* j_*)
/*Note -- these results are on the normalized scale -- adjusted later in program*/
matrix agg_results[2,1] = (e(tau_bc)-invnorm(.975)*e(se_tau_rb))*N_window/365
matrix agg_results[2,2] = e(tau_bc)*N_window/365
matrix agg_results[2,3] = (e(tau_bc)+invnorm(.975)*e(se_tau_rb))*N_window/365

/********************************************************************************
*** 4b. Parametric prediction                                                 ***
********************************************************************************/
/**NOTE: This uses a uniform kernel to get better estimates on the slopes to the left and right of the cutpoint.
This has minimal effect on the RD estimates. 
 ***/
capture drop s_hat*
capture drop diff
set matsize 1000
matrix param_results = J(500,1,.)
forval i=1(1)500 {
	preserve
	qui keep if abs(time_a)<45
	bsample _N
	qui reg sentence_normed_truncated announce##c.time_a i.judge i.statute_code if abs(time_a)<45
	qui predict s_hat1 if e(sample) & time_a>=0 & time_a<45
	/* Counterfactual predictions; assuming time trend of zero*/
	qui replace announce = 0
	qui replace time_a =0
	qui predict s_hat0 if e(sample) & time_a>=0 & time_a<45
	qui gen diff = max_possible*(s_hat1-s_hat0)
	qui summ diff
	scalar effect = r(mean)*r(N)/365
	disp in yellow effect
	matrix param_results[`i',1]=effect
	restore
}
capture drop param_results*
svmat param_results
centile param_results, centile(2.5 50 97.5)
matrix agg_results[3,1] = r(c_1)
matrix agg_results[3,2] = r(c_2)
matrix agg_results[3,3] = r(c_3)

/********************************************************************************
*** 4c. CIA-Based Estimators                                                  ***
********************************************************************************/
/*Determine sample based on positivity restrictions*/
preserve
keep if time_a > -40 & time_a < 45
egen positivity_judge = mean(announce), by(judge)
egen positivity_statute = mean(announce), by(statute_code)
summ positivity*
keep if positivity_judge > 0 & positivity_judge < 1 & positivity_statute > 0 & positivity_statute<1
count
/*Test Conditional independence assumptions */
qui reg sentence_normed_truncated i.judge i.statute_code numberofcounts time_a if announce == 0 
lincom time_a
qui reg sentence_normed_truncated  i.judge i.statute_code numberofcounts time_a if announce == 1 
lincom time_a
restore

/********************************************************************************
*** 4c(1). Linear Reweighting Estimator										  ***
********************************************************************************/
capture program drop linreweight
program  linreweight, rclass
	version 14.2
	syntax varlist(numeric fv), t(varname numeric) m(varname numeric)
	tokenize `varlist'
	qui reg `1' `2' `3' if `t' == 0
	tempvar s_hat0
	qui _predict `s_hat0' if `t' == 1, xb
	qui reg `1' `2' `3' if `t' == 1
	qui tempvar s_hat1
	qui _predict `s_hat1' if  `t' == 1, xb
	qui tempvar diff
	qui gen `diff' = (`s_hat1'-`s_hat0')*`m'
	qui summ `diff'
	return scalar d = r(mean)*r(N)/365
	scalar d1 = r(mean)*r(N)/365
	disp in yellow d1
end


matrix lrwe = J(500,1,.)
forval i=1(1)500 {
	preserve
	qui keep if time_a>-40 & time_a<45
	qui egen positivity_judge = mean(announce), by(judge)
	qui egen positivity_statute = mean(announce), by(statute_code)
	qui keep if positivity_judge > 0 & positivity_judge < 1 & positivity_statute > 0 & positivity_statute<1
	bsample _N
	gen one = 1
	linreweight sentence_normed_truncated i.judge i.statute_code numberofcounts, t(announce) m(max_possible)
	matrix lrwe[`i',1]=r(d)
	restore
	}	
capture drop lrwe
svmat lrwe
centile lrwe, centile(2.5 50 97.5)
matrix agg_results[4,1] = r(c_1)
matrix agg_results[4,2] = r(c_2)
matrix agg_results[4,3] = r(c_3)


/********************************************************************************
*** 4d. Propensity Score Estimator                                            ***
********************************************************************************/

matrix pse_results = J(500,1,.)
forval i=1(1)500 {
	preserve
	qui keep if time_a>-40 & time_a<45
	qui egen positivity_judge = mean(announce), by(judge)
	qui egen positivity_statute = mean(announce), by(statute_code)
	qui keep if positivity_judge > 0 & positivity_judge < 1 & positivity_statute > 0 & positivity_statute<1
	qui logit announce i.judge i.statute_code numberofcounts
	predict lambda, p
	bsample _N
	qui count if announce == 1
	local n_t1 = r(N)
	qui summ announce if e(sample)
	local tmean = r(mean)
	qui gen y_cf = sentence_normed_truncated*max_possible*(announce-lambda)/(`tmean'*(1-lambda)) 
	qui summ y_cf
	matrix pse_results[`i',1]=r(mean)*`n_t1'/365
	drop y_cf
	restore
	}
capture drop pse_results
svmat pse_results
centile pse_results, centile(2.5 50 97.5)
matrix agg_results[5,1] = r(c_1)
matrix agg_results[5,2] = r(c_2)
matrix agg_results[5,3] = r(c_3)


preserve
qui keep if time_a>-40 & time_a<45
tab announce
qui egen positivity_judge = mean(announce), by(judge)
qui egen positivity_statute = mean(announce), by(statute_code)
qui keep if positivity_judge > 0 & positivity_judge < 1 & positivity_statute > 0 & positivity_statute<1
tab announce
restore
	
/********************************************************************************
*** 4e. Graph Aggregate Results                                              ***
********************************************************************************/
capture drop agg_results*	
svmat agg_results
capture drop n
gen n = _n in 1/5
replace n = 0.5 in 6
replace n = 5.5 in 7
twoway  (rspike agg_results1 agg_results3 n , lcolor(gs11) lwidth(medthick)) (scatter agg_results2 n, msymbol(circle_hollow) mcolor(gs6) ), ///
	xtitle("") legend(off) xlabel(1 `" "(1a)" "Unadjusted LATE" "Extrapolation" "' 2 `" "(1b)" "Adjusted LATE" "Extrapolation" "' /// 
	3 `" "(2)" "Fully" "Parametric" "' 4 `" "(3)" "Linear" "Reweighting" "' 5 `" "(4)" "Propensity" "Score" "', ///
	labsize(medsmall)) ytitle(Additional Years Incarceration)	

graph export "aggregate_effects.pdf", replace	
l agg_results* in 1/6

	  
/********************************************************************************
*** 5. Additional Analyses	                                                  ***
********************************************************************************/

/********************************************************************************
*** 5a. Hetereogenous Effects by Leniency - Petition Only (Table C.1 in SI)  ***
********************************************************************************/

*residualize normalized sentence before the petition announcement
reghdfe sentence_normed_truncated  numberofcounts if time_a<0, absorb(statute_code judge_code, savefe)

*generate indicator whether judge is above median severity
sum  __hdfe2__, d //the coefficients associated with the second fixed effect i.e. judges
gen x = __hdfe2__>`r(p50)'
replace x = . if mi(__hdfe2__)
bys judge: egen above_median = max(x) 

/********************************************************************************
*** 5a(1). Above Median Leniency                                              ***
********************************************************************************/

matrix rd_results = J(6,6,.)
rdrobust sentence_normed_truncated time_a if abs(time_a)<180 & above_median==1, kernel(tri) vce(cluster cluster_judge_statute)  all
matrix b_temp = e(beta_p_l)
matrix v_temp = e(V_rb_l)
matrix rd_results[1,1]=  e(tau_bc)
matrix rd_results[2,1]= e(se_tau_rb)
matrix rd_results[3,1]=b_temp[1,1]
matrix rd_results[4,1]=sqrt(v_temp[1,1])		
matrix rd_results[5,1]=e(b_l)
matrix rd_results[6,1]=e(N_h_l)+e(N_h_r)
rdrobust sentence_normed_truncated time_a if abs(time_a)<180 & above_median==1,  kernel(tri) vce(cluster cluster_judge_statute) 	covs(numberofcounts s_* j_*) all
matrix b_temp = e(beta_p_l)
matrix v_temp = e(V_rb_l)
matrix rd_results[1,2]=  e(tau_bc)
matrix rd_results[2,2]= e(se_tau_rb)
matrix rd_results[3,2]=b_temp[1,1]
matrix rd_results[4,2]=sqrt(v_temp[1,1])		
matrix rd_results[5,2]=e(b_l)
matrix rd_results[6,2]=e(N_h_l)+e(N_h_r)

/********************************************************************************
*** 5b(1). Below Median Leniency                                              ***
********************************************************************************/

rdrobust sentence_normed_truncated time_a if abs(time_a)<180 & above_median==0, kernel(tri) vce(cluster cluster_judge_statute)  all
matrix b_temp = e(beta_p_l)
matrix v_temp = e(V_rb_l)
matrix rd_results[1,3]=  e(tau_bc)
matrix rd_results[2,3]= e(se_tau_rb)
matrix rd_results[3,3]=b_temp[1,1]
matrix rd_results[4,3]=sqrt(v_temp[1,1])		
matrix rd_results[5,3]=e(b_l)
matrix rd_results[6,3]=e(N_h_l)+e(N_h_r)
rdrobust sentence_normed_truncated time_a if abs(time_a)<180 & above_median==0,  kernel(tri) vce(cluster cluster_judge_statute) covs(numberofcounts s_* j_*) all
matrix b_temp = e(beta_p_l)
matrix v_temp = e(V_rb_l)
matrix rd_results[1,4]=  e(tau_bc)
matrix rd_results[2,4]= e(se_tau_rb)
matrix rd_results[3,4]=b_temp[1,1]
matrix rd_results[4,4]=sqrt(v_temp[1,1])		
matrix rd_results[5,4]=e(b_l)
matrix rd_results[6,4]=e(N_h_l)+e(N_h_r)


/********************************************************************************
*** 5b. Effect on Appellate Judges (Table C.2 in SI)                          ***
********************************************************************************/

use "recall_data_appellate.dta", clear
destring _all, replace
egen panel_code = group(panel)
qui tab panel_code, gen(j_)

egen appeal_division_code = group(appeal_division)
qui tab appeal_division, gen(d_)

capture drop time_a 
capture drop time_r
gen time_a = date(edate,"YMD")-date("06/06/2016","MDY")
gen time_r = date(edate,"YMD")-date("06/05/2018","MDY")


capture drop rd_results*
matrix rd_results = J(6,4,.)

/********************************************************************************
*** 5b(1). Petition Announcement                                              ***
********************************************************************************/
rdrobust prodefendant time_a if abs(time_a)<180, kernel(tri) vce(cluster panel_code) all
matrix b_temp = e(beta_p_l)
matrix v_temp = e(V_rb_l)
matrix rd_results[1,1]=  e(tau_bc)
matrix rd_results[2,1]= e(se_tau_rb)
matrix rd_results[3,1]=b_temp[1,1]
matrix rd_results[4,1]=sqrt(v_temp[1,1])		
matrix rd_results[5,1]=e(b_l)
matrix rd_results[6,1]=e(N_h_l)+e(N_h_r)
rdrobust prodefendant time_a if abs(time_a)<180, kernel(tri) vce(cluster panel_code) covs(j_*) all
matrix b_temp = e(beta_p_l)
matrix v_temp = e(V_rb_l)
matrix rd_results[1,2]=  e(tau_bc)
matrix rd_results[2,2]= e(se_tau_rb)
matrix rd_results[3,2]=b_temp[1,1]
matrix rd_results[4,2]=sqrt(v_temp[1,1])		
matrix rd_results[5,2]=e(b_l)
matrix rd_results[6,2]=e(N_h_l)+e(N_h_r)

/********************************************************************************
*** 5b(1). Recall Election	                                              ***
********************************************************************************/

rdrobust prodefendant time_r if abs(time_r)<180, kernel(tri) vce(cluster panel_code) all
matrix b_temp = e(beta_p_l)
matrix v_temp = e(V_rb_l)
matrix rd_results[1,3]=  e(tau_bc)
matrix rd_results[2,3]= e(se_tau_rb)
matrix rd_results[3,3]=b_temp[1,1]
matrix rd_results[4,3]=sqrt(v_temp[1,1])		
matrix rd_results[5,3]=e(b_l)
matrix rd_results[6,3]=e(N_h_l)+e(N_h_r)
rdrobust prodefendant time_r if abs(time_r)<180, kernel(tri) vce(cluster panel_code) covs(j_*) all
matrix b_temp = e(beta_p_l)
matrix v_temp = e(V_rb_l)
matrix rd_results[1,4]=  e(tau_bc)
matrix rd_results[2,4]= e(se_tau_rb)
matrix rd_results[3,4]=b_temp[1,1]
matrix rd_results[4,4]=sqrt(v_temp[1,1])		
matrix rd_results[5,4]=e(b_l)
matrix rd_results[6,4]=e(N_h_l)+e(N_h_r)
svmat rd_results


	  
/********************************************************************************
*** 6. Washington State Robustness Check  (Figure B.3 in SI)			      ***
********************************************************************************/
	  
clear
use "recall_data_washington.dta", 

**generate calendar
gen time_a = edate -  mdy(6,6,2016)
gen announce = time_a>=0
egen cluster_judge_statute =group(judge Statute_RCW)
set matsize 10000
encode Statute_RCW, gen(charge_code)
encode judge, gen(judge_code)
recode TOTALCONFINEMENTDAYS .=0
local vrange = "0.2(0.1)0.6"
local fitmodel = "qfit"
local fitmodel2 = "lpoly"
local ar = 0.67
/********************************************************************************
*** 6a. Graphical Analysis         										      ***
********************************************************************************/
/********************************************************************************
*** 6a(1). Unadjusted             										      ***
********************************************************************************/
preserve
rdbwselect sentence_normed_truncated time_a if abs(time_a)<180, kernel(tri) vce(cluster cluster_judge_statute)	
scalar bw = e(h_mserd)
keep if abs(time_a)<=bw
fastxtile n_all = time_a, nq(50)
gen n = .
replace n = n_all
collapse (mean) sentence_normed_truncated (median) time_a, by(n)
gen weight = max(0,(bw-abs(time_a))/bw)
twoway 	   (`fitmodel2' sentence_normed_truncated time_a if time_a <= 0, lcolor(gs8) lpattern(solid)) ///
	       (`fitmodel2' sentence_normed_truncated time_a if time_a >=0, lcolor(gs8) lpattern(solid)) ///
		   (`fitmodel' sentence_normed_truncated time_a if time_a <= 0 [aweight=weight], lcolor(maroon) lpattern(dash)) ///
	       (`fitmodel' sentence_normed_truncated time_a if time_a >=0 [aweight=weight], lcolor(maroon) lpattern(dash)) ///
		   (scatter sentence_normed_truncated time_a, mcolor(black) msymbol(circle_hollow)), ///
			name(graph_announce_unadjusted_wa, replace) ///
			title("Petition Announcement, Unadjusted")  xtitle("Date") ///
			xlabel( 0 "Jun 6" -45 "Apr 22" 45 "Jul 21",angle(45) nogrid) ///
			ytitle("Normalized Sentence") ylabel(`vrange', format(%3.1f) nogrid) ///
			xline(0, lcolor(black) lpattern(dash)) legend(off) scheme(plotplain) aspectratio(`ar') ///
			graphregion(color(white)) bgcolor(white)
graph export "WA_announce_normal_no_control.pdf", replace
restore

/********************************************************************************
*** 6a(2). Unadjusted             										      ***
********************************************************************************/
preserve	
rdbwselect sentence_normed_truncated time_a if abs(time_a)<180, kernel(tri) vce(cluster cluster_judge_statute)	
scalar bw = e(h_mserd)
qui reg sentence_normed_truncated i.judge_c i.charge_code  if abs(time_a)<=bw
qui predict ry if e(sample), resid
qui summ sentence_normed_truncated if e(sample), meanonly
qui replace ry = r(mean) + ry
qui reg time_a i.judge_c i.charge_code if abs(time_a)<=bw
qui predict rx if e(sample), resid
qui summ time_a if e(sample), meanonly
qui replace rx = r(mean)+rx 
fastxtile n_all = rx, nq(50)
gen n = .
replace n = n_all
collapse (mean) ry (median) rx time_a, by(n)
gen weight = max(0,(bw-abs(rx))/bw)
sum ry rx
twoway 	   (`fitmodel2' ry rx if rx <= 0, lcolor(gs8) lpattern(solid)) ///
	       (`fitmodel2' ry rx if rx  >=0 , lcolor(gs8) lpattern(solid)) ///
	       (`fitmodel' ry rx if rx <= 0 [aweight=weight], lcolor(maroon) lpattern(dash)) ///
	       (`fitmodel' ry rx if rx  >=0 [aweight=weight], lcolor(maroon) lpattern(dash)) ///
		   (scatter ry rx, mcolor(black) msymbol(circle_hollow)) ///
		    ,  name(graph_announce_adjusted_wa, replace) ///
			title("Petition Announcement, Adjusted")  xtitle("Date") ///
			xlabel( 0 "Jun 6" -45 "Apr 22" 45 "Jul 21",angle(45) nogrid) ///
			ytitle("Normalized Sentence") ylabel(`vrange', format(%3.1f) nogrid) ///
			xline(0, lcolor(black) lpattern(dash)) legend(off) scheme(plotplain) aspectratio(`ar') ///
			graphregion(color(white)) bgcolor(white)
			graph export "WA_announce_normal_control.pdf", replace
restore

gr combine graph_announce_unadjusted_wa graph_announce_adjusted_wa, ///
		   imargin(0 0 0 0) graphregion(margin(l=22 r=22)) name(all_rd_plots, replace)

graph export "combined_rd_plots_wa.pdf", replace



/********************************************************************************
*** 7. Descriptive Statistics												  ***
********************************************************************************/
gen male = arrest_gender_def=="Male"
replace male =. if arrest_gender_def=="NA"|mi(arrest_gender_def)

capture drop results*
matrix results = J(13,6,.)

local varlist = "sentence_days sentence_normed_to_statute sentence_normed_truncated sex_crime violent_nonsex_667 nonviolent_crime_667 Black Hispanic White male arrest_age_def"
forvalues i = 1/11 {
	local var : word `i' of `varlist'
	qui sum `var'
	matrix results[`i',1]=`r(mean)'
	matrix results[`i',2]=`r(sd)'
	qui sum `var' if abs(time_a)<45
	matrix results[`i',3]=`r(mean)'
	matrix results[`i',4]=`r(sd)'
	qui sum `var' if abs(time_r)<45
	matrix results[`i',5]=`r(mean)'
	matrix results[`i',6]=`r(sd)'
	}

*number of cases
gen nvals=1
sum nvals
matrix results[12,1]=`r(sum)'
matrix results[12,2]=`r(sum)'
sum nvals if abs(time_a)<45
matrix results[12,3]=`r(sum)'
matrix results[12,4]=`r(sum)'
sum nvals if abs(time_r)<45
matrix results[12,5]=`r(sum)'
matrix results[12,6]=`r(sum)'

*number of defendants
drop nvals
by defendant, sort: gen nvals = _n == 1 
replace nvals =. if mi(defendant)
sum nvals
matrix results[13,1]=`r(sum)'
matrix results[13,2]=`r(sum)'
sum nvals if abs(time_a)<45
matrix results[13,3]=`r(sum)'
matrix results[13,4]=`r(sum)'
sum nvals if abs(time_r)<45
matrix results[13,5]=`r(sum)'
matrix results[13,6]=`r(sum)'